using namespace std;
// PM has added this comment as a test of svn
#include<iostream> 
#include<fstream>
#include<cstdlib>
#include<math.h>

//#include "randomc.h"
#include"flysim.h"


cFly::cFly(cFlyCategory * flyCategory):cFlySim()
{
  iTotalOthers = gciTotalNumberOfFlies+1;
  this->flyCategory = flyCategory;
  //if (dPreComputedPowers == NULL){
   // dPreComputedPowers = new float[gciTotalNumberOfFlies+1];
 // }
}

cFly::~cFly()
{
}

//// More efficient way of calculating rates that we often recalculate
//float cFly::dRateMM(float dAttraction, int iHowManyMales)
//{
//	int i;
//	float dTol=1e-6;
//	static float dLastAttraction=-9;
//    //static float dPreComputedPowers[gciTotalNumberOfFlies+1]; // PPP: We need to put back this array, and prerably make it a single array that all flies share
//	if ((dAttraction-dLastAttraction>dTol)||(dLastAttraction-dAttraction>dTol))
//	{
//	  // "this is the first visit for a new value of the Male attraction to Males parameter, so reset all stored values\n";
//		for (i=0;i<gciTotalNumberOfFlies+1;i++)
//			dPreComputedPowers[i]=-9;
//		dLastAttraction=dAttraction;
//	}
//	
//	// Now try to lookup the value
//    //PPP: Once the array dPreComputedPowers[] is returned, we need to uncomment the line below
////	if (dPreComputedPowers[iHowManyMales]<-8) // need to calculate it for the first time
//        {
//		dPreComputedPowers[iHowManyMales]=pow(exp(dAttraction), log((float)(iHowManyMales))); //BRF_08_10
//        }
//        if (dPreComputedPowers[iHowManyMales] == 0){
//          cerr<<"DEBUG:"<<iHowManyMales<<","<<dAttraction<<endl;
//          cerr<<"DEBUG2: "<<exp(dAttraction)<<","<<log((float)(iHowManyMales))<<endl;
//          throw "MM Divide by zero";
//        }
//	return 1/dPreComputedPowers[iHowManyMales];
//}
//
//float cFly::dRateMF(float dAttraction, int iHowManyFemales)
//{
//	int i;
//	float dTol=1e-6;
//	static float dLastAttraction=-9;
//    //static float dPreComputedPowers[gciTotalNumberOfFlies+1]; // PPP: We need to put back this array
//	if ((dAttraction-dLastAttraction>dTol)||(dLastAttraction-dAttraction>dTol))
//	{
//		// this is the first visit for a new value of the Male attraction to Females parameter, so reset all stored values
//		for (i=0;i<gciTotalNumberOfFlies+1;i++)
//			dPreComputedPowers[i]=-9;
//		dLastAttraction=dAttraction;
//	}
//	
//	// Now try to lookup the value
//    //PPP: Once the array dPreComputedPowers[] is returned, we need to uncomment the line below
////	if (dPreComputedPowers[iHowManyFemales]<-8) // need to calculate it for the first time
//        {
//		dPreComputedPowers[iHowManyFemales]=pow(exp(dAttraction), log((float)(iHowManyFemales+1))); //BRF_08_10
//        }
//        if (dPreComputedPowers[iHowManyFemales] == 0){
//          cerr<<"DEBUG:"<<iHowManyFemales<<","<<dAttraction<<endl;
//          cerr<<"DEBUG2: "<<exp(dAttraction)<<","<<log((float)(iHowManyFemales+1))<<endl;
//          throw "MF Divide by zero";
//        }
//	return 1/dPreComputedPowers[iHowManyFemales];
//}
//
//float cFly::dRateFM(float dAttraction, int iHowManyMales)
//{
//	int i;
//	float dTol=1e-6;
//	static float dLastAttraction=-9;
//    //static float dPreComputedPowers[gciTotalNumberOfFlies+1]; // PPP: We need to put back this array
//	if ((dAttraction-dLastAttraction>dTol)||(dLastAttraction-dAttraction>dTol))
//	{
//		// this is the first visit for a new value of the Female attraction to Males parameter, so reset all stored values
//		for (i=0;i<gciTotalNumberOfFlies+1;i++)
//			dPreComputedPowers[i]=-9;
//		dLastAttraction=dAttraction;
//	}
//	
//	// Now try to lookup the value
//    //PPP: Once the array dPreComputedPowers[] is returned, we need to uncomment the line below
////	if (dPreComputedPowers[iHowManyMales]<-8) // need to calculate it for the first time
//		dPreComputedPowers[iHowManyMales]=pow(exp(dAttraction), log((float)(iHowManyMales+1))); //BRF_08_10
//        if (dPreComputedPowers[iHowManyMales] == 0){
//          cerr<<"DEBUG:"<<iHowManyMales<<","<<dAttraction<<endl;
//          cerr<<"DEBUG2: "<<exp(dAttraction)<<","<<log((float)(iHowManyMales+1))<<endl;
//          throw "FM Divide by zero";
//        }
//	return 1/dPreComputedPowers[iHowManyMales];
//}
//
//float cFly::dRateFF(float dAttraction, int iHowManyFemales)
//{
//	int i;
//	float dTol=1e-6;
//	static float dLastAttraction=-9;
//    //static float dPreComputedPowers[gciTotalNumberOfFlies+1]; // PPP: We need to put back this array
//	if ((dAttraction-dLastAttraction>dTol)||(dLastAttraction-dAttraction>dTol))
//	{
//		// this is the first visit for a new value of the Female attraction to Females parameter, so reset all stored values
//		for (i=0;i<gciTotalNumberOfFlies+1;i++)
//			dPreComputedPowers[i]=-9;
//		dLastAttraction=dAttraction;
//	}
//	
//	// Now try to lookup the value
//    //PPP: Once the array dPreComputedPowers[] is returned, we need to uncomment the line below
////	if (dPreComputedPowers[iHowManyFemales]<-8) // need to calculate it for the first time
//		dPreComputedPowers[iHowManyFemales]=pow(exp(dAttraction), log((float)(iHowManyFemales)));
//        if (dPreComputedPowers[iHowManyFemales] == 0){
//          cerr<<"DEBUG:"<<iHowManyFemales<<","<<dAttraction<<endl;
//          cerr<<"DEBUG2: "<<exp(dAttraction)<<","<<log((float)(iHowManyFemales))<<endl;
//          throw "FF Divide by zero";
//        }
//	return 1/dPreComputedPowers[iHowManyFemales];
//}
//

void cFly::dCalcRateOfLeaving()
{
  //cerr<<"cFly::dCalcRateOfLeaving\n";
  // Calculates the rate at which a fly wants to leave its present home	
	if (WhereAmI->iIndex!=0) 
	{
          //" I am currently on a petri dish, and I'm not sure I like the look of that fly over there - he's been staring at me"<<endl;
		switch (cSex)
		{
				//float dRate1=0;
				//float dRate2=0;
			case 'M': {
// this should not be zero or else we get overflows
assert(iTotalOthers+WhereAmI->iHowManyMalesAreHere>0);
				dLeaveRate=expRateLeaveMale*flyCategory->current_cache->dPreComputedPowers[cFlyCategory::INTERACTION_MM*iTotalOthers+WhereAmI->iHowManyMalesAreHere]*flyCategory->current_cache->dPreComputedPowers[cFlyCategory::INTERACTION_MF*iTotalOthers+WhereAmI->iHowManyFemalesAreHere];
				//dRate=dRate*dRateMM(MaleAttractionToMales,WhereAmI->iHowManyMalesAreHere);
				//cerr<<"I am fly "<<iID<<" I am on patch:"<<WhereAmI->iIndex<<". The number of males here is:"<<WhereAmI->iHowManyMalesAreHere<<"\n";
				//dRate=dRate*dRateMF(MaleAttractionToFemales,WhereAmI->iHowManyFemalesAreHere);
		//dRate=dRate*(1/(pow(MaleAttractionToMales, log((float)WhereAmI->iHowManyMalesAreHere))))*(1/(pow(MaleAttractionToFemales, log((float)WhereAmI->iHowManyFemalesAreHere + 1)))); 
			}
				break;
			case 'F': {
// this should not be zero or else we get overflows
assert(iTotalOthers+WhereAmI->iHowManyFemalesAreHere>0);
				dLeaveRate=expRateLeaveFemale*flyCategory->current_cache->dPreComputedPowers[cFlyCategory::INTERACTION_FM*iTotalOthers+WhereAmI->iHowManyMalesAreHere]*flyCategory->current_cache->dPreComputedPowers[cFlyCategory::INTERACTION_FF*iTotalOthers+WhereAmI->iHowManyFemalesAreHere];
				//dRate=dRate*dRateFM(FemaleAttractionToMales,WhereAmI->iHowManyMalesAreHere);
				//dRate=dRate*dRateFF(FemaleAttractionToFemales,WhereAmI->iHowManyFemalesAreHere);
				//cerr<<"I am fly "<<iID<<" I am on patch:"<<WhereAmI->iIndex<<". The number of females here is:"<<WhereAmI->iHowManyFemalesAreHere<<"\n";
		//dRate=dRate*(1/(pow(FemaleAttractionToMales, log((float)WhereAmI->iHowManyMalesAreHere + 1))))*(1/(pow(FemaleAttractionToFemales, log((float)WhereAmI->iHowManyFemalesAreHere))));
			}
				break;
				
				
			default: // should never get to this point
				throw"\nReached default state in switch statement in float cFly::dCalcRateOfLeaving(). Exit.";
		}
	}
	else 
	{
             //I am not on a petri dish, so I want to land on one\n";
		switch (cSex)
		{
			case 'M':
				dLeaveRate=expRateJoinMale;
				break;
			case 'F':
				dLeaveRate=expRateJoinFemale;
				break;
			default:
				// should never get to this point
				throw"\nReached default state in switch statement in float cFly::dCalcRateOfLeaving(). Exit.";
		}
	}
        //cerr<<"Fly leave rate: "<<dLeaveRate<<endl;
	//dLeaveRate=dRate;
}


void cFly::dOldCalcRateOfLeaving()
{
	// Calculates the rate at which a fly wants to leave its present home	
	float dRate;
	if (WhereAmI->iIndex==0) // I am not on a petri dish, so I want to land on one
	{
		switch (cSex)
		{
			case 'M':
				dRate=expRateJoinMale;
				break;
			case 'F':
				dRate=expRateJoinFemale;
				break;
			default:
				// should never get to this point
				throw"\nReached default state in switch statement in float cFly::dCalcRateOfLeaving(). Exit.";
		}
	}
	else // I am currently on a petri dish, and I'm not sure I like the look of that fly over there - he's been staring at me
	{
		switch (cSex)
		{
				//float dRate1=0;
				//float dRate2=0;
			case 'M': {
				dRate=expRateLeaveMale;
				//dRate1= log((float)WhereAmI->iHowManyMalesAreHere);
				//dRate2= log((float)WhereAmI->iHowManyFemalesAreHere + 1);
				dRate=dRate*(1/(pow(exp(flyCategory->current_cache->dAttractionConstant[cFlyCategory::INTERACTION_MM]), log((float)WhereAmI->iHowManyMalesAreHere))))*(1/(pow(exp(flyCategory->current_cache->dAttractionConstant[cFlyCategory::INTERACTION_MF]), log((float)WhereAmI->iHowManyFemalesAreHere + 1)))); 
			}
				break;
			case 'F': {
				dRate=RateLeaveFemale;
				//dRate1= log((float)WhereAmI->iHowManyMalesAreHere + 1);
				//dRate2= log((float)WhereAmI->iHowManyFemalesAreHere);
				dRate=dRate*(1/(pow(exp(flyCategory->current_cache->dAttractionConstant[cFlyCategory::INTERACTION_FM]), log((float)WhereAmI->iHowManyMalesAreHere + 1))))*(1/(pow(exp(flyCategory->current_cache->dAttractionConstant[cFlyCategory::INTERACTION_FF]), log((float)WhereAmI->iHowManyFemalesAreHere))));
			}
				break;
				
				
			default: // should never get to this point
				throw"\nReached default state in switch statement in float cFly::dCalcRateOfLeaving(). Exit.";
		}
	}
	dLeaveRate=dRate;
}

void cFly::Print(bool bNewLine)
{
	if (bNewLine)
		ofile<<endl;
	ofile<<"ID: "<<iID<<",  Sex: "<<cSex<<",  Genotype: "<<iGenotype;
	if (WhereAmI==NULL)
		ofile<<",  Petri Dish: I have no location.";
	else
		ofile<<",  Petri Dish: "<<WhereAmI->iIndex;
}

int cFly::LandOnAPetriDish()
{
	// for now we just pick one at random - it could of course depend on presence of other flies, food, etc.
	int iDish;
//	float p;
	//p=gsl_rng_uniform_pos(rg);
	//iDish=(int) 1+floor(p*gciNumberOfPetriDishes);
	// the 0-th petri dish is off site so that is why we add one.
	iDish=(int)1+gsl_rng_uniform_int(rg,gciNumberOfPetriDishes);
        //if (iDish>gciNumberOfPetriDishes) iDish=gciNumberOfPetriDishes;
	return iDish;
}

int cFly::LeaveThisPetriDish(int iDishBeingLeft)
{
	// A fly decides to leave its current petri dish - for now we assume it goes to the external environment
	int iNewHome;
	iNewHome=0;
	return iNewHome;
}

float cFly::getMaleAttractionToMales(){
  return flyCategory->current_cache->dAttractionConstant[cFlyCategory::INTERACTION_MM];
}

float cFly::getMaleAttractionToFemales(){
  return flyCategory->current_cache->dAttractionConstant[cFlyCategory::INTERACTION_MF];
} 

float cFly::getFemaleAttractionToMales(){
  return flyCategory->current_cache->dAttractionConstant[cFlyCategory::INTERACTION_FM];
}

float cFly::getFemaleAttractionToFemales(){
  return flyCategory->current_cache->dAttractionConstant[cFlyCategory::INTERACTION_FF];
}

void cFlyCategory::setMaleAttractionToMales(float val){
  float eval = exp(val);
  if (val!=current_cache->dAttractionConstant[INTERACTION_MM]){
    //cerr<<"Recomputing for MM\n";
    for(int others=0;others<iTotalOthers;++others){
      current_cache->dPreComputedPowers[INTERACTION_MM*iTotalOthers+others]=
      1/pow(eval, dPowers[others]);
    }

    //cerr<<"Recomputed\n";
  }
  current_cache->dAttractionConstant[INTERACTION_MM] = val;
} 

void cFlyCategory::setMaleAttractionToFemales(float val){
  float eval = exp(val);
  if (val!=current_cache->dAttractionConstant[INTERACTION_MF]){
    //cerr<<"Recomputing for FM\n";
    for(int others=0;others<iTotalOthers;++others){
      current_cache->dPreComputedPowers[INTERACTION_MF*iTotalOthers+others]=
      1/pow(eval, dPowers[others+1]);
    }
    //cerr<<"Recomputed\n";

  }
  current_cache->dAttractionConstant[INTERACTION_MF] = val;
}
void cFlyCategory::setFemaleAttractionToMales(float val){
  float eval = exp(val);
  if (val!=current_cache->dAttractionConstant[INTERACTION_FM]){
    //cerr<<"Recomputing for FM\n";
    for(int others=0;others<iTotalOthers;++others){
      current_cache->dPreComputedPowers[INTERACTION_FM*iTotalOthers+others]=
      1/pow(eval, dPowers[others+1]);
    }

    //cerr<<"Recomputed\n";
  }
  current_cache->dAttractionConstant[INTERACTION_FM] = val;
}

void cFlyCategory::setFemaleAttractionToFemales(float val){
  float eval = exp(val);
  if (val!=current_cache->dAttractionConstant[INTERACTION_FF]){
    //cerr<<"Recomputing for FF\n";
    for(int others=0;others<iTotalOthers;++others){
      current_cache->dPreComputedPowers[INTERACTION_FF*iTotalOthers+others]=
      1/pow(eval, dPowers[others]);
    }
    //cerr<<"Recomputed\n";
  }
  current_cache->dAttractionConstant[INTERACTION_FF] = val;
}

cache_container_t::cache_container_t(int iTotalOthers){
  dPreComputedPowers = new float[4*iTotalOthers];
  dAttractionConstant = new float[4];
  dPreComputedInteractions = new float[4*iTotalOthers*2];  
}

cache_container_t::~cache_container_t(){
  delete[] dPreComputedPowers;
  delete[] dAttractionConstant;
  delete[] dPreComputedInteractions;
}

void cFlyCategory::copyCache(cache_container_t * in, cache_container_t * out){
//cerr<<"copying\n";
  for(int i=0;i<4;++i){
    for(int j=0;j<iTotalOthers;++j){
      out->dPreComputedPowers[i*iTotalOthers+j] = 
      in->dPreComputedPowers[i*iTotalOthers+j];
    }
    out->dAttractionConstant[i] = in->dAttractionConstant[i];
  }
  for(int i=0;i<iTotalOthers*2;++i){
     //out.dPreComputedInteractions[] = in.dPreComputedInteractions[i];
  }
//cerr<<"copied\n";
}

void cFlyCategory::commitParameters(){
  //cerr<<"Committing\n";
  copyCache(current_cache,previous_cache);
  //cerr<<"Committed\n";
}

void cFlyCategory::rollbackParameters(){
  //cerr<<"Rolling back\n";
  copyCache(previous_cache,current_cache);
  //cerr<<"Rolled back\n";
}

cFlyCategory::cFlyCategory(){
  iTotalOthers = gciTotalNumberOfFlies+1;
//  cerr<<"totalothers for "<<iTotalOthers<<"\n";
  dPowers = new float[iTotalOthers+1];
  for(int others=0;others<=iTotalOthers;++others){
    dPowers[others] = log(static_cast<float>(others));
  }
  current_cache = new cache_container_t(iTotalOthers);
  previous_cache = new cache_container_t(iTotalOthers);

}
cFlyCategory::~cFlyCategory(){
  delete current_cache;
  delete previous_cache;
  delete[]dPowers;
}
